{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PDE-FIND for identifying Navier-Stokes\n",
    "\n",
    "Samuel Rudy, 2016\n",
    "\n",
    "This notebook demonstrates PDE-FIND for the vorticity equation given a simulation of fluid flowing around a cylinder.\n",
    "\\begin{align*}\n",
    "\\omega_t &= \\frac{1}{Re}\\nabla^2 \\omega - (V \\cdot \\nabla)\\omega\\\\\n",
    "V &= (v,u)\\\\\n",
    "Re &= 100\n",
    "\\end{align*}\n",
    "The x and y components of the velocity field are given as forcing terms to the PDE.  That is, they appear in $\\Theta$, but are not differentiated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "pylab.rcParams['figure.figsize'] = (12, 8)\n",
    "import numpy as np\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from PDE_FIND import *\n",
    "import scipy.io as sio\n",
    "import itertools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load in the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Load data\n",
    "data = sio.loadmat('../CYLINDER/DATA/ALL.mat')\n",
    "steps = 151\n",
    "n = 449\n",
    "m = 199\n",
    "W = data['VORTALL'].reshape(n,m,steps)   # vorticity\n",
    "U = data['UALL'].reshape(n,m,steps)      # x-component of velocity\n",
    "V = data['VALL'].reshape(n,m,steps)      # y-component of velocity\n",
    "\n",
    "dt = 0.2\n",
    "dx = 0.02\n",
    "dy = 0.02\n",
    "\n",
    "# Cut out the portion of the data before the cylinder\n",
    "xmin = 100\n",
    "xmax = 425\n",
    "ymin = 15\n",
    "ymax = 185\n",
    "\n",
    "W = W[xmin:xmax,ymin:ymax,:]\n",
    "U = U[xmin:xmax,ymin:ymax,:]\n",
    "V = V[xmin:xmax,ymin:ymax,:]\n",
    "n,m,steps = W.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we take the SVD of the data and reconstruct either with a reduced basis or everything.  It isn't necesarry but is interesting to show that we can still derive the correct PDE with the first 50 modes (maybe less). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "W = W.reshape(n*m,steps)\n",
    "U = U.reshape(n*m,steps)\n",
    "V = V.reshape(n*m,steps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "uw,sigmaw,vw = svd(W, full_matrices=False); vw = vw.T\n",
    "uu,sigmau,vu = svd(U, full_matrices=False); vu = vu.T\n",
    "uv,sigmav,vv = svd(V, full_matrices=False); vv = vv.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fe6e57b5f50>]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs8AAAHhCAYAAACRPFNsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl0VfW5//H3NwlhEAQZZJ5HRXACBGUIg+BQpa1D1Vo7\n3bbaW+1o7Ww623pre9tabX+1Xsdaa2vRqgioRwUcQBRUxjLLDCKKzGH//thBQ0gw5JycfZK8X2ud\nBWfvnb0fs+xaH7999vMNURQhSZIk6YPlJV2AJEmSVFsYniVJkqQqMjxLkiRJVWR4liRJkqrI8CxJ\nkiRVkeFZkiRJqiLDsyRJklRFhmdJkiSpirISnkMITUIIs0II52TjeZIkSVJNyNbK83XA37L0LEmS\nJKlGHHF4DiHcFkLYEEKYV+74WSGEhSGExSGE68ocHwfMBzYBIe2KJUmSpISEKIqO7AdCGA5sB+6M\nomhg6bE8YDEwFlgLzAIuiaJoYQjhJ0AToD+wI4qij2SwfkmSJClrCo70B6Iomh5C6Fru8BBgSRRF\nKwFCCPcBE4GFURR9r/TYFcDmiu4ZQjiyBC9JkiRVUxRF1e6GyFTPc0dgdZnvb5Qee08URXdGUfRo\nZTeIoshPhj7XX3994jXUlY+/S3+fufzx9+nvM1c//i79febyJ12OqpMkSZKqKFPheQ3Qpcz3TqXH\nJEmSpDqjuuE5cPDkjFlArxBC1xBCIXAJ8FC6xal6ioqKki6hzvB3mVn+PjPL32dm+fvMHH+XmeXv\nM7dUZ9rGvUAR0ArYAFwfRdHtIYSzgd8QB/Lboii64QjuGV1//fUUFRX5L4gkSZIyLpVKkUql+OEP\nf0iUxguDRxyea0IIIcqFOiRJklS3hRDSCs++MChJkiRVkeFZkiRJqiLDsyRJklRFhmdJkiSpinIm\nPBcXF5NKpZIuQ5IkSXVQKpWiuLg47fs4bUOSJEn1htM2JEmSpCwxPEuSJElVZHiWJEmSqsjwLEmS\nJFWR4VmSJEmqIsOzJEmSVEWGZ0mSJKmKciY8u0mKJEmSaoqbpEiSJElHyE1SJEmSpCwxPEuSJElV\nZHiWJEmSqsjwLEmSJFWR4VmSJEmqIsOzJEmSVEWGZ0mSJKmKDM+SJElSFRmeJUmSpCoyPEuSJElV\nZHiWJEmSqsjwLEmSJFWR4VmSJEmqopwJz9dffz2pVCrpMiRJklQHpVIpiouL075PiKIo/WrSLSKE\naMuOLbRs3DLpUiRJklSHhRCIoihU9+dzZuV53Tvrki5BkiRJOqycCc/rt69PugRJkiTpsHImPK/b\n7sqzJEmSclvOhOflmw3PkiRJym05E55v+9s69u5NugpJkiSpcjkTnvc0XMdVV0EODP+QJEmSKpQz\n4bnXSet5+WX46U+TrkSSJEmqWEHSBRywccc6nvo3DBsG+/fDGWdA9+7QuTM0aJB0dZIkSVIWVp5D\nCP1CCLeEEO4PIVxZ2XXrtq+jfXuYPBmWL4cf/xjGjIGmTWH8+JquUpIkSfpgWdthMIQQgDuiKLqi\ngnNRwx835M3r3qRJgyYHnduzB9q0gWXLoFWrrJQqSZKkOirrOwyGEG4LIWwIIcwrd/ysEMLCEMLi\nEMJ15c6dB/wbeLSy+7Zr2q7CjVIKC+Hkk+Gll460UkmSJCmzqtO2cTswoeyBEEIe8PvS4/2BS0MI\n/Q6cj6Lo4SiKzgUur+ym7Zu1r3SL7lNPNTxLkiQpeUf8wmAURdNDCF3LHR4CLImiaCVACOE+YCKw\nMIQwCvgo0BB4pLL7tm/avtJdBk89Ff75zyOtVJIkScqsTE3b6AisLvP9DeJATRRFTwNPf9ANVk9a\nzW1Tb+O1Tq9RVFREUVHRe+dOPRW+970MVSpJkqR6I5VKkUqlMna/ar0wWLry/HAURQNLv18ATIii\n6POl3y8HhkRRdE0V7xf9KPUjdu7byc/G/uyQ8/v3Q4sWsGIFtGx5xOVKkiRJQAIvDFZiDdClzPdO\npceqrH2z9hW+MAiQlwcnnQRz5lS/QEmSJCld1Q3PofRzwCygVwihawihELgEeOhIbni4nmfwpUFJ\nkiQlrzqj6u4FZgJ9QgirQgifjqKoBLgamAK8DtwXRdGCI7nv4aZtgOFZkiRJyavOtI3LKjn+GPBY\ndQu557f3sPKtlZWeP/VUuP766t5dkiRJ9VmmXhzM2g6Dhy0ihGhvyV4a/7QxO7+7k4K8QzN9SUn8\n0uCqVXDMMQkUKUmSpFovV14YTFtBXgGtGrdi47sbKzyfnx+/NPjyy1kuTJIkSSqVM+EZ7HuWJElS\nbsut8PwBEzdOOcXwLEmSpOTkVHhu17SdK8+SJEnKWZnanjttxR070mtIe7Z9qg2cWvE1/frBunWw\nbRs0b57d+iRJklR71blpG9HkySz6zffo/Ow8mnzzu/CDH1R47RlnwE9+AqNHZ7lISZIk1Xp1ZtoG\nEybw2q++xXd/cAbcc0+ll516qtt0S5IkKRm5E56Jp228ePR2WL0aduyo8Br7niVJkpSUnArP7Zq2\n441dG6BPH3j99QqvMTxLkiQpKTnzwiDEo+rWb19PNHAU4dVXYfDgQ67p1w82bYJzzoEePaBnT+jV\nCyZMgMLCBIqWJElSvZFT4blxg8Y0LmjMzuN60WTevAqvKSiAuXPjz9Kl8eeuu+Dhh+FPf8pywZIk\nSapXcio8Q9z3vLFNO7pNTVV6TefO8eeAd96BIUPgttvgs5+t+RolSZJUP+VUzzPErRururaAefOg\nimP0mjWDBx+Eb30LZs+u4QIlSZJUb+XMynNxcTFFRUW0a9qOlQ13QV5evCNKhw5V+vl+/eDWW+HC\nC+MA3bp1DRcsSZKkWqPubZJSWsfXH/867Zq249rrH4dvfAPOOuuI7nXddfEc6MmTIT+/JqqVJElS\nbVV3Nkkp1b5Ze9ZtXwcDB8atG0fopz+N/+zZE8aPhyuvhBtvhCVLMlyoJEmS6p2cads4oH3T9jyy\n5BG29Z5I85lHPtC5oAAeewyWLYsncSxbBjNmwLRp8PjjNVCwJEmS6o2ca9vY9O4mvjz5y6xNPcz/\n+9d+/vHX73Nx/4vpcUyPat9/48Z435U334xbqSVJklQ/1bm2jTZHteHeC+5lys/foMemfbyxaRkj\nbh+R1j2PPTZ+gXDBggwVKUmSpHop58LzAYVNm5PfvQe/6/kl3t79Nlt3bk3rfsOGwXPPZag4SZIk\n1Us5G54BGDiQ8Oqr9G7ZmyVvpvfGn+FZkiRJ6cr58My8efRp1YclWwzPkiRJSlatCM+9W/Zm8ZbF\nad1qwABYvTp+aVCSJEmqjtwOzwMGwKuvxivPabZtFBTAoEHwwgsZqk2SJEn1Tm6H565d4Z13OC7v\n2LRXnsHWDUmSJKUnt8NzCDBgAH3W7mbJm0tIdya14VmSJEnpyO3wDDBwIEcvWkGDvAZsfHdjWrca\nOhRefBFKSjJUmyRJkuqVWhGemTeP3q3SH1fXpg20bQvz52eoNkmSJNUrOROei4uLSaVSh54YOBDm\nzqVPqz72PUuSJKlaUqkUxcXFad8npNtHnAkhhKjSOnbsgA4d+NWdV7GpScQN425I61m33hpP3Lj9\n9rRuI0mSpFoohEAURaG6P58zK8+VatIEzj6botmb027bAFeeJUmSVH25H54BLr2UPlPnZKRt44QT\nYO1a2LIlA3VJkiSpXqkd4XnCBJouXs6u5f9hf7Q/rVvl58PgwfD88xmqTZIkSfVG7QjPDRsSPvwR\nLl/QgDfefiPt29m6IUmSpOooSLqAKrv0Uj72ub+zZMsSujTvktathg2Da66BnTuhaVM46ijo0QMu\nvDBDtUqSJKlOqh0rzwCjR9P+rX1seHl62rcaNw6+/W3o0AHy8mDzZrjySlicfku1JEmS6rDcH1VX\nxksXDWdj4/2cfefMjNfwhS9A797wjW9k/NaSJEnKEXV/VF0Z7370PPo/8SrUQOA//3x46KGM31aS\nJEl1SI2H5xDCxBDCn0IIfw0hnJnOvVqP/RB5u3bBvHmZKu89Y8fC3LlxC4ckSZJUkRoPz1EUTYqi\n6PPAVcDF6dyrZ8te/LV/xP6/3puZ4spo1CgO0I8+mvFbS5IkqY444vAcQrgthLAhhDCv3PGzQggL\nQwiLQwjXVfCj3wNurm6hAA0LGvLU0LaU3HN3PKh569Z0bneI88+HSZMyektJkiTVIdVZeb4dmFD2\nQAghD/h96fH+wKUhhH5lzt8APBpF0Stp1ApAyYD+rB5/Glx9NXTpAm3bwujRMH9+urfm3HNh2jTY\ntSvtW0mSJKkOOuLwHEXRdKD8ku8QYEkURSujKNoL3AdMBAghXA2MBS4MIXw+zXrp3aoPD39uFMya\nBW+/DXPmwIQJ8JnPQElJWvdu0wYGDoSnnkq3SkmSJNVFmdokpSOwusz3N4gDNVEU/Q743QfdoLi4\n+L2/FxUVUVRUVOF1fVr1YdHmRfGXEKBjR/jmN+Nm5Vtvhf/+7+r9E5Q6MHXj7LPTuo0kSZJyQCqV\nIpVKZex+1ZrzHELoCjwcRdHA0u8XABNKXwwkhHA5MCSKomuqeL8qzXkGeGzJY/z6+V8z5RNTDj6x\nYAGMHAmvvBIH6mpatAjGjIE33oizuSRJkuqOXJnzvAYou2d2p9JjGdenVR8Wb6lgK8DjjoMvfjHu\nhU5D377QrBm89FJat5EkSVIdVN3wHEo/B8wCeoUQuoYQCoFLgBrZcqRri66s376e7z/5fW5/+Xae\nXvE0a94uzenf/ja8/jr8619pPcMNUyRJklSR6oyquxeYCfQJIawKIXw6iqIS4GpgCvA6cF8URQuO\n5L7FxcVV6kcpyCvggYsfoCCvgKdWPMV3n/wux918HE+veDoe1vzHP8arz2+/faT/aO8xPEuSJNUt\nqVTqoHfsqqtaPc+ZdiQ9zxX59rRv07CgIcVFxfGBL30J7rwTTjkFBg+OP2eeCcccU6X7lZRA+/bx\njzVpAo0bx39+85vQo0e1y5QkSVLC0u15rhPhedLCSdz60q089vHH3j+4ZQvMnh1/Hn8cOnWCe6u+\nM+GSJbBwIezcGX8eeQS6doUbb6x2mZIkSUqY4RlYv309/f/Qn83XbiZUNCLj5ZfhE5+A116r9jOe\new6uvBLmzq32LSRJkpSwXJm2kah2TdvRtLApS95cUvEFffvC0qWwb1+1nzF4MKxaBevXV/sWkiRJ\nquXqRHgGOK3jabzwxgsVn2zSJG5iXras2vcvKIh3AZ82rdq3kCRJUi1XZ8Lz0E5DeWFNJeEZ4Pjj\nYf78tJ4xfjxMmfLB10mSJKluypnwXNVRdZU5reNphw/Pxx0X70KYhgPhOQfaxCVJknQEHFVXzs69\nO2n1y1Zs+eYWGjdofOgFt98OTz4Jd92V1nN69YJ//hMGDkzrNpIkSUqALwyWatygMce1OY6X179c\n8QUZWHmGeFz01Klp30aSJEm1UJ0JzwBDOw6t/KXB446LBzfv35/WM+x7liRJqr/qVHg+rdNpPL/m\n+YpPNm8ef1avTusZo0fDzJnxximSJEmqX+pWeD7cuDrISOtGixZxv/P06WndRpIkSbVQnQrPvVv1\nZtvubWzYvqHiCzIwrg7i1g37niVJkuqfOhWe80IeQzoOqXxkXQZfGrTvWZIkqf7JmfCc7pznA4Z2\nHMrzb1TS95yhlechQ2DFCrfqliRJqi2c81yJR5c8yq+e+xVPXPHEoSc3boR+/WDLFgjVHu8HwEc/\nCoMHxy0cBQXQoAF06gRHH53WbSVJklSD0p3zXOfC8+Ydm+n52568+c03yc/LP/hkFEGbNvDaa9Cu\nXVrPeegh+MlPYN8+2LsXdu2CZs1gzpy0bitJkqQa5CYp5bRu0po2TdqwcPPCQ0+GkLG+5/PPhxdf\njMPyq6/CokWweXP8d0mSJNVNBUkXUBOGdhrKVY9cRd9WfWla2JRmDZsxvud4hncZ/n54Hj06o8/M\ny4PLLoN77oEbbsjorSVJkpQj6lzbBsCqbauYsWoG7+x5h+17tvP6xtdZuW0l066YBr/5DfznP/D7\n32fseQe8+iqce278MmFenVvTlyRJqv3SbduokyvPXZp3ocuALu993/juRvr+vi9RFBGOOw4efrhG\nnjtgABxzDDz7LIwaVSOPkCRJUoLqxfrosUcdy9ENj2bp1qUZG1dXmY9/PG7dkCRJUt1TL8IzwKAO\ng5i1ZlY8T277dti6tUaec+ml8I9/wO7dNXJ7SZIkJaj+hOf2g5i9dnZGJ25UpHNnGDgQHn20Rm4v\nSZKkBNWb8Dy442Bmr5sdf6nB8Axx68bdd9fY7SVJkpSQehOeT21/KnPWzaFkf0nc91yD4fnCC2Ha\nNHjrrRp7hCRJkhJQb8LzMY2Poe1RbVm0ZVG88jxjBrzwAqxcGW8PmEEtWsC4cXHvsyRJkuqOOjmq\nrjKDOw5m9trZHD/s7Hh77muugfXr40+XLvGg5kaNMvKsj38cfvzjePvuBg2goCAeY/ehD8Vt15Ik\nSap96lV4PvDS4BUnXgEPPvj+iSiCQYNg9mwYPjwjzzr3XJg1C156KQ7Q+/bBk09C06YZ39xQkiRJ\nWZIz4bm4uJiioiKKiopq7BmDOgzigQUPHHoihHhXk2eeyVh4btgQfv7zg4/deCPcd5/hWZIkKdtS\nqRSpVCrt+9TJ7bkr887ud2j3q3a8dd1bNMhvcPDJf/0Lbr0VJk+useevXAmnngpr10JhYY09RpIk\nSZVId3vuevPCIECzhs3o2rwr8zdVsMPg8OHw3HNxf0UN6doV+vaNJ3FIkiSp9qlX4RlKdxpcO+vQ\nE61bx7sPzp1bo8+/9NK4dUOSJEm1T70Lz4M7xBM3KjRiRNz3XIMuuggefhh27qzRx0iSJKkG1Lvw\nPKjDoMrD88iR8OyzNfr8tm3jwR5u3y1JklT71LvwfGK7E5m/aT679+0+9OSBlecafnnxkkvgr3+t\n0UdIkiSpBtS78NykQRP6tOrDvA3zDj3ZuTM0a1ajW3cDfPSjMHUqvP12jT5GkiRJGVbvwjMk37px\nzDHxWOlJk2r0MZIkScoww3N5I0fW+EuDELduOHVDkiSpdqnxTVJCCN2B7wJHR1F0cSXXZGWTlANe\n2/gaJ//xZAryCijML6Qwv5BjjzqWl7/wMoVLV8DYsbBqVbzzYA3Zvh06downbzRvHu9IWFgYz4LO\nz6+xx0qSJNVr6W6SkrUdBkMI9+dKeAbYt38fe0v2sqdkD3v372Xk7SO5+6N3c0q7k6F9e3j+eejW\nrUZr+MlP4taN3bthzx7YsgU+/3n46U9r9LGSJEn1VtZ3GAwh3BZC2BBCmFfu+FkhhIUhhMUhhOuq\nW1C2FOQV0LhBY5o3ak7rJq0Z0nEIs9bMilebs9S68b3vwaxZMG8eLFwY7zx4zz01PuxDkiRJ1VSd\nnufbgQllD4QQ8oDflx7vD1waQuhX7udqrgciAw7aPCULm6VUZODAuH1jVgUbIEqSJCl5Rxyeoyia\nDmwtd3gIsCSKopVRFO0F7gMmAoQQWoYQbgFOyuUV6cEdB7+/bXcWJm5UJIR4B8K//z3rj5YkSVIV\nFGToPh2B1WW+v0EcqImi6E3gqg+6QXFx8Xt/LyoqoqioKEOlVc2JbU9k8ZbF7Ni7gyYnnAAbN0Jx\nMXToAMceG39OPTVeGq5BF10EEyfCL39Zo+8rSpIk1QupVIpUKpWx+1XrhcEQQlfg4SiKBpZ+vwCY\nEEXR50u/Xw4MiaLomireL+svDFZk0J8G8duzf8vpnU+HyZNh+vQ4RG/cCHPmwHe+A1deWaM1RBH0\n6wd33QVDhtTooyRJkuqddF8YzNTK8xqgS5nvnUqP1SqDOwxm1ppZcXg+66z4c8Af/gAvvVTjNZRt\n3TA8S5Ik5ZbqbpISOPgFwFlArxBC1xBCIXAJ8FC6xWXbQX3P5Z18Mrz8clbqOBCec2AxXpIkSWVU\nZ1TdvcBMoE8IYVUI4dNRFJUAVwNTgNeB+6IoWpDZUmve4A6HCc8DB8L8+bB3b43XMXBgvGGKUzck\nSZJyyxG3bURRdFklxx8DHqtuIcXFxYm8KFjWcW2OY83ba3hr11u0aNTi4JNHHRVv/7dgQZxua1AI\ncPHFtm5IkiRlSqZeHMzaDoOHLSJHXhgEGP6X4fyw6IeM7TH20JOXXQYTJsAnP1njdcydG0/dWL7c\nqRuSJEmZkvUdBuu6gzZLKS+Lfc+2bkiSJOWeTE3bqDMGdxzMPxf8s+KTJ58M//53Vuo40LrxzW/C\nsGFxkG7QIA7V55+flRIkSZJUjivP5Rz2pcGTT4ZXXoH9+7NSy9e+Fgflo4+GvDzYsQMuvxzefTcr\nj5ckSVI5rjyX06tlL97e/TYb393IsUcde/DJVq2geXNYtgx69arxWlq2jAN0WS+8ANOmxf3QkiRJ\nyq6cWXkuLi7O6NaJ1RVCYFCHQcxak/y854qcdx48/HBij5ckSaqVUqkUxcXFad/HaRsV+M4T36Ew\nv5DiouJDT/7wh7B7N/zsZ1mvC2DpUjjjDFi7Nm7lkCRJUtU5baMGfGDfc4Irzz17xt0jTuGQJEnK\nPsNzBQ60bVS4Gp5weAZbNyRJkpJi20YFoiiiw00dGNl1JI0LGlOQV0BBXgFfGfoV+rXqC23awKuv\nQvv2idQ3YwZ88YvxRiqSJEmqOts2akAIgX9e/E/O63MeRd2KOK3jaWx8dyN3zr0zHsCc8Orz0KFx\nz/PKlYmVIEmSVC85qq4SwzoPY1jnYe9979aiGz965kfxlwPh+ZxzEqktPx/OPTdu3fjSlxIpQZIk\nqV5y5bmKhnUexpx1c9i5d2fiK89g37MkSVISDM9V1LSwKf3b9OfFNS/mRHgePx6eew7eeSfRMiRJ\nkuqVnAnPubJJyuGM7DqSZ1c9C717w8aNsG1bYrU0awbDhsGUKYmVIEmSVGu4SUoCJi2cxM2zbmbK\nJ6bA6afDz38Oo0YlVs/NN8OLL8IddyRWgiRJUq2S7rQNw/MR2LJjC93/tztvXvcmBVd/GXbuhHHj\n4jf48vLgxBOhT5+s1bNyJZx0EnziE1BQAA0awFFHwVe/Gq9MS5Ik6WCG5yw74Q8ncPvE2xm8al+8\n9FtSAvv3w7vvwpIlsGhRVuv5179g1SrYuzf+/Otf8NnPwuc+l9UyJEmSagXDc5Z98ZEv0qtlL742\n7GsHn9i/Hzp0gJkzoUePZIoD/vpXuPdeJ3FIkiRVxE1SsmxElxE8s/KZQ0/k5cFZZ8Hkydkvqoyz\nz4ann44XwiVJkpRZhucjNKLrCKavms7+aP+hJ3MgPLdoAUOGwNSpiZYhSZJUJxmej1CnozvRvFFz\nFmxacOjJM8+Ml313785+YWWcfz489FCiJUiSJNVJhudqGNFlRDzvubxWreD442H69OwXVcb558O/\n/x2/yyhJkqTMMTxXw8iuIyvue4a46fixx7JbUDndukH79vD884mWIUmSVOcYnqvhwEuDFU4IOeus\nxMMz2LohSZJUEwzP1dCrZS9KohJWvLXi0JODBsVbd69alfW6ypo4ESZNSrQESZKkOsfwXA0hBEZ0\nGcE/FvyDRZsXsfTNpax8ayXb92yPR9aNH5/41I1TToF33sn6ni2SJEl1mpukVNMD8x/gB0/9gJKo\nhH3797Fr3y66tejGjM/MgLvvhn/+M/4k6KqroGdP+MY3Ei1DkiQpZ7jDYI7YtW8XrX/ZmjVfW0Pz\nt3dDnz5x+0ZhYWI1PfYY/Oxn8GwFg0EkSZLqozqzw2BxcTGpVCrpMqqtUUEjTut0WjzC7thjoXdv\neO65RGsaPRrmzYNNmxItQ5IkKXGpVIri4uK07+PKcwb99JmfsmXnFm6acBN8//uwdy/ccEOiNV16\nKbzwAjRvDg0axJ8vfAGuuCLRsiRJkhKR7spzQSaLqe/GdB/DVY9cFX85+2y45BLYvx/y8+NPy5Zw\n9dVxgs2S226DpUvjHL93LyxcCDfeaHiWJEmqDleeM2hvyV5a39iapdcspXWjlvB//webN8db/e3b\nB3fdBTffHG/jnZCSEmjXDl56Cbp0SawMSZKkRLjynEMa5DdgeJfhpFakuPD4C+Eznzn4grw8ePTR\nRMNzfv77+7h84QuJlSFJklQr5cwLg3XFmG5jeHL5kxWfPOecODwn7Nxz4ZFHkq5CkiSp9jE8Z9iY\n7ocJzyedFO9c8p//ZLeocsaPh1QKdu1KtAxJkqRax/CcYSe2O5FNOzax5u01h54MIX6RMOHV55Yt\nYeBAeOaZRMuQJEmqdQzPGZYX8ijqVsRTK56q+IJzz008PEPcQWLrhiRJ0pExPNeAw/Y9jxsHM2bA\nu+9mt6hyciTDS5Ik1So1Hp5DCE1CCP8XQvhjCOGymn5eLhjTfQxPLH+CCsfvHX00DB4MT1YSrrNk\n4EDYsQOWLEm0DEmSpFolGyvPHwX+HkXRF4Dzs/C8xPVr3Y89JXtY/tbyii/IgakbIdi6IUmSdKSO\nODyHEG4LIWwIIcwrd/ysEMLCEMLiEMJ1ZU51AlaX/r0kjVprjRDC4aduHAjPCW8MkwMZXpIkqVap\nzsrz7cCEsgdCCHnA70uP9wcuDSH0Kz29mjhAA1R7N5fa5rB9z8cdF2+YMn9+dosqZ9w4eO452L49\n0TIkSZJqjSMOz1EUTQe2ljs8BFgSRdHKKIr2AvcBE0vPPQhcGEK4GXg4nWJrk7E9xvLgwgfp+/u+\nDLhlAKf+6VQm3D2Bt3a9lTM9E82awWmnwRNPJFqGJElSrZGp7bk78n5rBsAbxIGaKIp2AJ+p6IfK\nKi4ufu/vRUVFFBUVZai0ZHRr0Y3/XP0ftu/Zzp6SPezdv5evT/k6U5dO5aL+F8Xh+cYb4ZvfTLTO\nc86B226Lt+1u2BAKC6FTJ+jZM9GyJEmSMiKVSpFKpTJ2v1DhRIgP+qEQugIPR1E0sPT7BcCEKIo+\nX/r9cmBIFEXXVPF+UXXqqG1+8/xvWLBpAX8874/xqIt27eJxF23axG0cCVi3Dr761bh1Y/fu+PPK\nK7B6NTRvnkhJkiRJNSaEQBRF1W4lzlRiWwN0KfO9U+kxlXFmjzOZumxq/KVJE7joIujeHQoK4mXf\nFi3gL38WABs8AAAgAElEQVTJak3t28N998G//w1Tp8a7Dp5+OkybltUyJEmSaoXqhufAwS//zQJ6\nhRC6hhAKgUuAh9Itrq45vs3x7Nq3i6VvLo0P3HZbvAJdUgJvvw3/939w++2J1gjxDuKPPZZ0FZIk\nSbmnOqPq7gVmAn1CCKtCCJ+OoqgEuBqYArwO3BdF0YLMllr7hRAY12Mc05ZNK38iXnmeMAHmzoW3\n3kqmwFJnnw2TJyc+SU+SJCnnVGfaxmVRFHWIoqhhFEVdoii6vfT4Y1EU9Y2iqHcURTcc6X2Li4sz\n2sydqw5q3SivcWMYMSLun0hQ795xln/11UTLkCRJyphUKnXQgIrqqtYLg5lWX14YBFj3zjr6/6E/\nm67dRH5e/qEX3HwzzJ6dePvG1VfHUzeuu+6Dr5UkSaotcuWFQVVR+2bt6dCsA3PWzan4ggM9E/v3\nZ7ewCspw90FJkqSDGZ4TcNjWjR494hlxr7yS3aLKKSqCOXNg27ZEy5AkScophucEjOsxrvLwDDkx\n7qJJEzjjDEfWSZIklWV4TsCobqOYvXY27+55t+ILcqRnIgcyvCRJUk7JmfBcX6ZtADQtbMop7U/h\n2VXPVnzByJHxqIs338xuYeU4sk6SJNUVTtuo5X789I95a9db/GrCryq+4Lzz4OMfh0suyW5hZUQR\n9OoFDz4IAwcmVoYkSVLGOG2jljqz55lMW36YhuJzzkm8ZyIEWzckSZLKMjwnZFCHQazatoqZq2cy\nf9N8FmxawKLNi9i9b3d8gSPrJEmSco5tGwm6dsq1PPqfR4miiP3Rfrbu2spnT/4sPxv7s/iC44+H\nO++EQYMSq3HHDmjbFr74xXgCx1FHQdOmcUdJs2aJlSVJklQt6bZtGJ5zyPRV0/ny5C/z0udfig98\n/euweTNMnAgNGkBBAXTuDCeckNW6HnsM5s6Ng/SOHfD00/Cxj8E3vpHVMiRJktJmeK5D9pbspfWN\nrVl6zVJaN2kNixfD978Pe/bAvn3xny+9BBs3Ql5yHTeTJsW7iE+ZklgJkiRJ1WJ4rmM+dO+H+OSJ\nn+Si/hdVfEHfvvC3v8FJJ2W3sDLefhs6dowzfOPGiZUhSZJ0xJy2UceM6zGOacsOM4Vj9Gh46qns\nFVSBo4+Os/szzyRahiRJUtblTHiuT5ukHM64HuMOP8IuB8IzwPjxtm1IkqTaw01S6qgoiuhwUwdm\nfGYGPY7pcegFGzbErRubN8cvECbkxRfhs5+NN0KUJEmqLWzbqGNCCIdv3WjbNm44fvnl7BZWzqmn\nwtq18UeSJKm+MDznoHHdP6DvecyYxFs38vNh7FhbNyRJUv1ieM5BY3uM5cnlT7I/qmR3QfueJUmS\nEmF4zkGdju5Em6Pa8Mr6Vyq+YNQomDED9u7NbmHljB8PU6cmvoO4JElS1hiec9RhWzdatYIePWDW\nrOwWVU6XLtC6NbxSScaXJEmqawzPOerMnmfmfN8z2LohSZLqF8NzjhrVdRTPvfEcu/btqviCHOp7\nfvzxpKuQJEnKDsNzjmreqDkDjh3AjFUzKr5g5Eh44QXYvTu7hZUzahTMng3btydahiRJUlYkt8uG\nPtC4HuP4zQu/Yf6m+eTn5ZMX8uhxTA/G9xwPzZtDv37w/PNxgk1I06YweDBMngwXXACh2iPHJUmS\ncp8rzznskyd+kh4terBoyyJe2/gac9bN4aK/X8Sekj3xBTnS9/yZz8AXvgDNmsHAgfDhD8P99ydd\nlSRJUua5PXctM+hPg7hpwk2M7DoyXu79+c/h6aeTLguAbdtg+XKYMweuvx5WrXIlWpIk5ZZ0t+c2\nPNcy33niO+SHfH485sdxo3HbttCrFzRqFH+aN4dbb4UOHRKrMYqgZ0+YNAkGDEisDEmSpEOkG55z\npm2juLiYVCqVdBk578weZzJlWelsuKZNYdEiuPNO+O1v4Yc/jJd6//GPRGsMAc4+Gx57LNEyJEmS\n3pNKpSguLk77Pq481zK79+2mzY1tWPmVlRzT+JhDL/jb3+Cee+Chh7JfXBn//jf86lc50ZItSZL0\nnjqz8qyqaVjQkOFdhvPk8icrvmDs2LgHOuGtu0ePjkfYvf12omVIkiRllOG5Fjqzx5lMWVrJtn6t\nW8c90M8/n92iyjnqKDj9dJh2mE0SJUmSahvDcy00vud4pi6bWvkFZ54JUw9zPkvse5YkSXWN4bkW\nOr7N8ewu2c3SN5dWfMH48TkVnm1nlyRJdYXhuRYKITCux7jKWzfOOANeew22bs1uYeX06QOFhXEp\nkiRJdYHhuZYa3+MwrRsNG8YBOuFRFyHAOefAo48mWoYkSVLGGJ5rqXE9xvHUiqfYt39fxRfY9yxJ\nkpRxhudaqm3TtnRp3oVZa2ZVfEGOhOfRo+Gll+KtuyVJkmo7w3MtdtjWjQED4u27ly3LblHlNGni\nyDpJklR31Gh4DiF0DyH8OYRwf00+p746s+dh5j2HkDOrz+ecY+uGJEmqG7KyPXcI4f4oii4+zHm3\n566GnXt30u1/u7Fv/z4K8wtpkNeAxg0aM+mSSfRr3Q/uvDPepvuBBxKtc+lSGDwYfvc7uOyyONdL\nkiQlId3tuasUnkMItwEfAjZEUTSwzPGzgN8Qr2DfFkXRLyr5ecNzDdm1bxc79u5gT8ke9pTsoThV\nzHGtj+PaM66FtWvhhBNg0ybIz0+0zhdfhM99Dtq3h1tuge7dEy1HkiTVU+mG56q2bdwOTCj34Dzg\n96XH+wOXhhD6lZ77RAjhphBC+wOXV7dAHV6jgka0bNySdk3b0aV5Fyb2nci05aUNxh06xJ+774YX\nXojf3Js7F95+O+t1DhkCs2dDUVG8Cv3rX7t5iiRJqn0KqnJRFEXTQwhdyx0eAiyJomglQAjhPmAi\nsDCKoruAu0IILUMItwAnhRCuq2xlGqC4uPi9vxcVFVFUVHRE/yCKFXUr4vIHL2fXvl00KmgE//3f\n8Kc/wb598Wfr1vgNvrvvznptDRrAt74FF14Io0bFLdknnJD1MiRJUj2SSqVIpVIZu1+Ve55Lw/PD\nB9o2QggXABOiKPp86ffLgSFRFF1zxEXYtpFRw24bxs/G/IzR3UcfenL58jg8r12baPPxlVdC797w\n9a8nVoIkSaqHstW2oVpkXPdxTFtWyWy47t2hcWNYsCC7RZUzYQJMqWRQiCRJUq5KJzyvAbqU+d6p\n9JgSNq7HuPf7nisyZgw8+WT2CqqkhJkzYefORMuQJEk6IkcSngMHv/g3C+gVQugaQigELgEeymRx\nqp6hnYayYNMCtu7cWvEFORCemzeHE0+EZ59NtAxJkqQjUqXwHEK4F5gJ9AkhrAohfDqKohLgamAK\n8DpwXxRF1e4FKC4uzmgzd33WsKAhp3c+ndSKVMUXjB4NqRSUlGSzrENMmACPP55oCZIkqZ5IpVIH\nDaiorqxskvKBRfjCYMb9z8z/YfnW5dx87s0VX3D88fHEjVNOyW5hZbzwAnz2s/Daa4mVIEmS6hlf\nGFSFakPf86BBsG4drLFTXpIk1RKG5zpqYNuBvLnzTVZtW1XxBTkQnvPzYexYp25IkqTaw/BcR+WF\nPMZ2H8sTy56o+IKiIpg+HfbsyWpd5TmyTpIk1SaG5zrssK0bLVvGu5TMmpXdosoZPx6mTk383UVJ\nkqQqyZnw7LSNzBvXI94spdKXMXOgdaNzZzj2WHj55UTLkCRJdZzTNlQlvX7bizs+fAcD2g6gYX5D\nCvMLCQe25X7sMfjlL+GppxKt8atfhdat4bvfTbQMSZJUDzhtQ4f1qZM+xfn3nU/HmzrS7OfNyPtR\nHl97/GvxyeHD47aNhLf5Gz/evmdJklQ7uPJcz8zbMI8L7r+AJVcviQ+ccQb86Efx2IuE7NgBbdvC\nypVxK7YkSVJNceVZR+SEY09g265t74+wy4G+5yZN4LLLoGNH6NULzjsPrrsOVlUyZU+SJCkprjzX\nQ5c8cAln9TqLT530KZg5E0aOhIYNoaAg/rRvD6+8Ev89i/btg6VLYf58uPfeuIzf/jarJUiSpDrO\nlWcdsTHdx/Dk8tLV5tNPh3ffhY0b46XexYvj4y+9lPW6Cgqgb1/4yEfg2mvB4SuSJCnXGJ7roTHd\nx/DE8ifeH2HXsCEcdRQ0bw6tWsG4cfBEJZurZMkpp8RZftOmRMuQJEk6iOG5Hup5TE8K8gpYvGVx\nxReMHZt4eC4oiIeBPP10omVIkiQdJGfCs5ukZE8I4eDWjfJGjYIXXkh8hN3o0YmPoJYkSXWEm6Qo\nLXfNvYtJiybxwMUPVHzBsGHwk58kOsJuzhy4/PL4BUJJkqRM8IVBVcuY7mNIrUixP9pf8QU50Pd8\n4omwbh2sX59oGZIkSe8xPNdTHY/uSOsmrZm3YV7FF+RA33N+fjxFz24eSZKUKwzP9diY7mN4Ylkl\nAXnYsLhf4q23sltUOfY9S5KkXGJ4rsfGdh/LkysqeWmwYcM4QCc87sLwLEmSconhuR4r6lbE9FXT\n2Vuyt+ILcqB1Y8AAePNNWLMm0TIkSZIAw3O91qpJK3oc04NZa2dVfEEOhOe8vHhynqvPkiQpFxie\n67mx3cdWPu/55JPjcRfr1mW3qHJs3ZAkSbnCOc/13OT/TOaSBy6h+zHdaZDXgML8Qjoe3ZH7LriP\nEAJ89KNwwQXw8Y8nVuPrr8N558GyZYmVIEmS6oh05zwbnuu5KIpYsHkBu/btYm/JXvaU7OHj//w4\nUz4xhX6t+8HNN8NLL8Ff/pJgjdCuHbz4InTtmlgZkiSpDkg3PBdkshjVPiEEjm9z/EHHxvYYS2pF\nKg7PY8fCL34RJ9hQ7X/P0qwRiori1o1PfSqREiRJkgB7nlWBoq5FPLWitMm4b19o0gTatoXu3eH4\n42HIEHjllazWNHEiXHklHHtsvPp83HHwX/+V1RIkSZJs29ChVry1gtP+fBrrv74+7nvetSveLGXn\nTtixA269FRo3hl/+Mqt1vfNO/PidO2HbNjjjDNiwAY46KqtlSJKkWizdtg1XnnWIbi260aRBExZu\nXhgfaNQobjru3h3694eLLkpk/EWzZvECeLducOKJ8ef557NehiRJqsdyJjwXFxeTSqWSLkOliroV\nkVqRqvjkaafBggXx8m+CRo5MfANESZJUS6RSKYqLi9O+j20bqtCdc+/k34v/zf0X3V/xBWPHwle/\nCh/6UHYLK2PyZLjhBvC/uSRJUlXZtqEaMarrKFIrUlT6HzU5sHPJGWfA7NlxS7YkSVI2GJ5Voa4t\nutK0sCkLNi+o+IIcCM/NmsVTN2ZVsru4JElSphmeVanD9j0PHgxLlsDWrVmtqbxRo+x7liRJ2WN4\nVqWKupWZ91xeYSEMGwbPPJPdosoZOTLxEiRJUj1ieFaliroV8fSKp3O673nEiHhc3d69iZYhSZLq\nCcOzKtWleReaNWzG/E3zK74gB8LzMcdAjx7w0kuJliFJkuoJw7MOq6jrYfqeTz0Vli+HLVuyWlN5\ntm5IkqRsMTzrsA7b99ygQTwvLuE39nxpUJIkZUuNh+cQwsQQwp9CCH8NIZxZ089TZhV1K+LplU+z\nY+8O9kf7D70gB1o3RoyAGTOgpCTRMiRJUj1Q4+E5iqJJURR9HrgKuLimn6fM6ty8M71b9qbVL1tR\n8KMCGv2kES1/0ZKpS6fGF4wenfgWf8ceCx06wNy5iZYhSZLqgSqH5xDCbSGEDSGEeeWOnxVCWBhC\nWBxCuO4wt/gecHN1C1VyZn52Jju/u5OSH5Sw9bqtfG3Y13h48cPxyZNPhtWrYdOmRGu071mSJGXD\nkaw83w5MKHsghJAH/L70eH/g0hBCv9Jznwgh3BRC6BBCuAF4NIqiVzJUtxIQQqBxg8ac2eNMnl5Z\n2mRcUADDh8M//pHoPtn2PUuSpGwIlc7wrejiELoCD0dRNLD0+1Dg+iiKzi79/i0giqLoF2V+5mrg\nCmAW8EoURX+q4L7RkdShZO0t2UurX7Zi+ZeX06pJK5g0Cb75TVi5Mu6h6N0bPvlJuOKKrNW0di30\n7QvnngsNG8afo46C4mJo3jxrZUiSpBwXQiCKolDdny9I8/kdgdVlvr8BDCl7QRRFvwN+90E3Ki4u\nfu/vRUVFFBUVpVmaakqD/Aac0eUMnl31LB/u92GYODH+7NsXt3A8/zxcdx184hMQqv3v5hHp0AEe\nfDDuHtm9O/786U8wdSpceGFWSpAkSTkolUqRyuD7WemuPF8ATCh9IZAQwuXAkCiKrjmiIlx5rnVu\nmH4D67ev5zdn/ebQk1EEXbrAk0/Gq9AJ+eUv4Y034Le/TawESZKUY9JdeU532sYaoEuZ751Kj6mO\nOzDCrkIhwLhxMG1adosqx5cIJUlSph1peA6lnwNmAb1CCF1DCIXAJcBDmSpOuevU9qey9M2lvLnz\nzYovGDs28fB8yimwdCls3ZpoGZIkqQ45klF19wIzgT4hhFUhhE9HUVQCXA1MAV4H7ouiaEHNlKpc\n0iC/AcM6D+PZlc9WfMHYsfHmKQnuXFJYCKedFm+gIkmSlAlVDs9RFF0WRVGHKIoaRlHUJYqi20uP\nPxZFUd8oinpHUXRDdQspLi7OaDO3at6orqMqb91o3z5+i2/OnOwWVY6tG5IkCeIXB8sOqKiuI3ph\nsKb4wmDtNHP1TL706JeY84VKAvJXvgJt28K3v53dwspIpeBb34oHgEiSJCX9wqDqsUEdBrHkzSW8\nteutii/IgZcGTzsNXn0Vtm9PtAxJklRHGJ5VbYX5hQztNLTyvudRo+CFF2DHjuwWVkbjxvEO4q48\nS5KkTDA8Ky1FXQ8zsq5ZMzjppMTf2LPvWZIkZYrhWWkZ1W0UqRWpyi/IgdYNw7MkScqUnAnPTtuo\nnQZ3GMyiLYvYtmtbxRfkwLzn00+H2bPjLbslSVL95LQN5Yyxd47lq0O/yof6fOjQk3v2QOvWsGxZ\n/GdCBg2C3/wGhg9PrARJkpQDnLahxF18/MVccP8FdP/f7oy7cxxfePgLTFtWutpcWAgjRsQbpiRo\n1ChbNyRJUvpceVZG7CnZw8q3VrJ061JmrZnFPa/ew8IvLYxP/vrXsHAh/PGPidU3aRLccgtMnpxY\nCZIkKQeku/JseFbG7Y/20+bGNsy7ch4dj+4Ir78Oo0fD3/4W/5mALVugR4/4z4KCREqQJEk5wLYN\n5Zy8kMeY7mN4YvkT8YH+/eH//T/45Cfhiitg48as19SqFXTpAn/+MzzxBEyfDrNmwbZK3nOUJEmq\niOFZNWJs97Hv9z0DTJwI8+dDmzZwwgnwl79kvaZrronbN37603jL7ksvhW98I+tlSJKkWsy2DdWI\nJVuWUHRHEW989Q1CKPf/jLzyCowfHy//9umTTIHAvHlw4YWweHFiJUiSpCyzbUM5qVfLXhTkFbBo\ny6JDT550Epx9Njz5ZPYLK+OEE2DzZli/PtEyJElSLZIz4dlNUuqWEMKhrRtljRmT+Pi6vLx47vOz\nzyZahiRJygI3SVHOu2fePTyw4AEe/NiDh55ctSreuWT9+jjFJuR//gdWroTf/S6xEiRJUhbZtqGc\nNab7GFIrUpTsLzn0ZJcu0Lx5PMYuQSNHunmKJEmqOsOzakz7Zu3p2KwjL617qeILRo9OvO/55JPj\nncPffDPRMiRJUi1heFaNGtt9LE8se6Lik2PGJB6eGzSAoUNhxoxEy5AkSbWE4Vk1amyPse9vllLe\n6NFxz0RJBW0dWWTrhiRJqirDs2rUqK6jeP6N59m5d+ehJ9u2hY4d4eWXs19YGYZnSZJUVYZn1ajm\njZozoO0AZq6eWfEFOdC6MWRI/N7i9u2JliFJkmoBw7Nq3Njuh2ndyIHw3Lhx/OLg888nWoYkSaoF\nDM+qceN6jOOhRQ/x1q63Dj05alT8tt6ePdkvrAxbNyRJUlUYnlXjhncZzoguIzju5uO4e97dHLQh\nzjHHQJ8+8OKLyRWI4VmSJFWNOwwqa1544wWufORKjml0DH849w/0a90vPnHttXD00fD97ydW2zvv\nQPv2sGULNGyYWBmSJKmGucOgao3TOp3GrM/NYmLfiZx+2+ls2L4hPpEDfc/NmkG/fjBrVqJlSJKk\nHGd4VlYV5BXw5aFfZlS3UUxbNi0+OHx4nFp3VjDOLots3ZAkSR/Etg0l4pZZt/D8mue548N3xAcm\nTIAlS+Ckk2DAgPhz7rnxKIwsefhhuOyyePR048bxp3dvuOOOrJUgSZJqWLptGwWZLCYdxcXFFBUV\nUVRUlHQpyoLxPcfz42d+TBRFhBDgkUdg8WJ49dX489OfwmuvQXFx1mr60Ifiec87dsSL4Dt3woc/\nDCtXQteuWStDkiTVgFQqRSqVSvs+rjwrMT1/25NJl0zihGNPOPTk5Mnwi1/AU09lv7AyLr4YzjkH\nPvWpRMuQJEkZ4guDqrXO7HEmU5ZOqfjk6afHfdC7d2e3qHJGj4YM/EeqJEmqIwzPSsz4nuMrD89H\nHw19+8Ls2dktqpzRo+PFb/+PEUmSBIZnJWhM9zHMXD2TXft2VXxBDoy/6Ns33vxw+fJEy5AkSTnC\n8KzEtGjUgv7H9mfGqhkVX5AD4TmE91efJUmSDM9K1Pgeh2ndGD4cnnsOSkqyW1Q5RUX2PUuSpJjh\nWYka33M8U5dNrfhkmzbx0OW5c7NbVDn2PUuSpAMMz0rUkI5DWLZ1GRvf3VjxBSNGJN660atX/Od/\n/pNoGZIkKQcYnpWoBvkNKOpW9P5W3eXZ9yxJknJIjYbnEEK/EMItIYT7QwhX1uSzVHsdtnVjxAh4\n9tnEeyYMz5IkCWo4PEdRtDCKoquAjwGn1+SzVHsdmPdc4S6TnTtDs2awcGH2CyvjwGYp9j1LklS/\nVSk8hxBuCyFsCCHMK3f8rBDCwhDC4hDCdZX87HnAv4FH0y9XdVHPY3pSmF/IVyZ/hdvm3MYzK59h\n7Ttr3w/TOdC60a0bFBbCokWJliFJkhJW1ZXn24EJZQ+EEPKA35ce7w9cGkLoV3ruEyGEm0II7aMo\nejiKonOByzNYt+qQEAJ/v+jvHHvUsTy76lm+88R36P+H/vz+xd/HF+TAS4P2PUuSJIBQ4f9VXtGF\nIXQFHo6iaGDp96HA9VEUnV36/VtAFEXRL8r8zCjgo0BDYG4URbdUcu+oqnWofpi6dCo/SP2A5z77\nHCxZAmPGwKpVcYpNyB13wCOPwP33J1aCJElKUwiBKIqqHSgK0nh2R2B1me9vAEPKXhBF0dPA01W5\nWXFx8Xt/LyoqoqioKI3SVNuN6jaKhZsXsn77etr16gX79sHKlXH/REJGj4Zrr4VXXolbOBo2hCZN\noH37xEqSJEkfIJVKkcrgbmfprDxfAEyIoujzpd8vB4ZEUXTNERfhyrMqcOk/LmVMtzF87tTPwcc+\nBueeC1dckWhNF10EixfDnj2wezds2gR33gkf+UiiZUmSpCpKd+U5nWkba4AuZb53Kj0mZcTEvhOZ\ntGhS/GX0aLjmGujbF046CYYNg69/Pes1/f3v8YaHCxbAsmVw001w331ZL0OSJCXkSFaeuxGvPA8o\n/Z4PLALGAuuAF4FLoyhacMRFuPKsCmzbtY3Ov+7M2q+vpWlBE3jjDdi5M/68/Taccw5s2RL3TyRk\nw4Y4z69fD40aJVaGJEmqoqysPIcQ7gVmAn1CCKtCCJ+OoqgEuBqYArwO3Fed4HxAcXFxRvtRVPs1\nb9ScoZ2GMmXpFMjLgy5d3l95Hjky/vvs2YnW2LYtDBgATzyRaBmSJOkDpFKpg96xq64qrzzXJFee\nVZmbX7yZF9e+yB0fvuPQk1/5Svy23nUVjhjPmptugvnz4c9/TrQMSZJUBUn2PEs17vy+5/PI4kfY\nt3/foSeHD4fp07NfVDkf+Qg89BCUlCRdiSRJqmmGZ+W0zs0707VFV2asmnHoyTPOgBkzYP/+7BdW\nRvfu0KFDXIokSarbDM/KeQdN3SirfXto2TLumUjYRz4CDz6YdBWSJKmmGZ6V8w6E5wr74keMyJnW\njQcfBFv3JUmq2wzPynkD2w6kZH8Jr296/dCTOdL3PGAA5OfHuw9KkqS6K53tuTOquLjYbblVoRAC\nE/tO5HtPfo8zOp9Bfl4++SGf49scz5nDh8OPf5x0iYTw/urzyScnXY0kSSovU9t0O6pOtcLqbau5\nZfYt7C3ZS0lUwr79+/jra3/l2U89Q7/jR8KcOdC5c6I1zpgBV10F8+YlWoYkSTqMdEfVGZ5Va/36\nuV/z+NLHeewfjQkXXwyXXppoPfv3vz91o2fPREuRJEmVcM6z6q0vDfkSq7atYn7fljnR95yXBxdf\nHLdhH398vBHikCHwox8lXZkkScoUV55Vq01dOpVb/vBpHpjakrwc6JfYuxdWroQ9e+LPjh3wsY/B\nXXeB7fySJCXPtg3VexfeM5F7PjeZhms3QIsWSZdziIcegm98A+bOhcaNk65GkqT6zbYN1Xu/POfX\nvNhhP5unPZR0KRU6//y4hSMHhoJIkqQ0ufKsOuGpK0ay5e31tLjpDzTIa0BBXgENCxrSv01/GjdI\nfrl3/XoYOBCmToUTT0y6GkmS6q90V56d86w6Yegl17LiG5/lS9N/zr79+9hbspcde3ewbOsyRnQd\nwTm9zuHs3mfTvUV3Qqj2/16qrV07uOEG+K//gueeg4Kc+V+eJEn1g3OepbLeeQfat49TaYMGUFgI\nTZqwbdL9PB4t4dEljzL5P5PZ8O4GAPJCHoEQ/xniP8seq+x4utf+Z0keTY8KtGiRR37pNfl5eTQ9\nKo+8vFCl+5Y9XqVra/if84jryaFrs/XvgaT/397dx1hWlwcc/z6z010QCxUMUEBYQFYF1wpNqCDo\nKFAQFDSmFlSsEEOsLWhrWoGGMJi0wYaWmPQlUZQoQYloC6g0UKSjUVxZslreFkHedoEChbYi6q7s\n3Kd/nDPs3evcmXPf5tx75/tJTu45v/Nyn31yd+7zO/ec35GGhzcMSnO2boVf/rIY5uKFF+DjHy+G\nuKqC2HcAAA6YSURBVPjwh3fYLDNpZINGNkia5sv2+dr6sW0jG2x+LPnrv2mwZWsy2yjatmxJnnq6\nwSFrGqxdmxx6WIOX7togIiEaTEwkGQ0iGjCRxWuU6ycaJOU2NIiJBIr1lPu/OE+DjAaUsVWJNwex\nLUvwHkO0LVC56B6Gwr91uyDmXVdpu063b7PdoI/R6fGXy3tL48riWWrnyivhppvgmmvqjmRRzz0H\n69YVD1hZt65YbjRgdraY5uZbX7tZF1GMSb1ixfbX5vl2r3Pzc1Prcru2xaa5eAa1vu73gKKTMzFR\ndGaaOzjxYlvRGcqWjlHEr29L2QEiGmWnKnfoIM3NJ9uP1WiaL9qTzKLj1Vz4zxX/rR2EJH+tw9jJ\ndoM+xnzrOjrGEvybul3Xj1x0+t7Nnb5BFfqD7kAN9L3HtLO0IlawYmIFkxOTrIjytcPl3XbajZUr\nVtb2XVqVxbPUzsMPw9FHwxNPFJWNgO4L8dlZyNxx3dzUutzaNrdfu6nX9cP2Hs25Wmi/+dp7aWtt\nXygGKP5bNBf9rR2AbtZ5jPE4RkQSc79oNXfQyo7cXEesdV3RsVv6DkgdnYyhfW+6O9ZsY5bZnGW2\nMcu2xjZms3ztYPnaP7iWqdVTdX7FVTI2NwxKfbd6dXH98wMPwJo1dUczNLafEdVylrm9oG4trKsW\n4L1uN0rH2LZtOOJYumMEjUbQaEx0dIzMhTtlw9xhqOMYEwusW+gXwKrr5psmJ9uvW2ianCxuJ/KG\nd4tnjbMIePObYWbG4llq0VzgSP3SrlM2+p2JwWw3O7u9Yzbf1MtlegtN27Ytvs18+2zdWvzdWLWq\nKKRXrdpx/oor4Kij6v4UDp7Fs8bb1BTceiucc07dkUjS2LNTNv7miuhf/ap4bZ4/4IC6o1saXvOs\n8faTnxQF9ObNXvcsSZJ6vubZvqHG28EHF68PPlhvHJIkaSxYPGu8zV33/O1v1x2JJEkaAxbPGn8W\nz5IkqU8snjX+pqaKETe8rl6SJPXI4lnj75BDisd1P/JI3ZFIkqQRZ/Gs8ed1z5IkqU+Gpnienp5m\nZmam7jA0ruYu3ZAkScvSzMwM09PTPR/HcZ61PGzcCCefDA8/XHckkiSpRo7zLFXx6lfDL34Bjz5a\ndySSJGmEWTxreYiAE06A174W9tqreIbomjVwxRV1RyZJkkaIl21o+Zidheeeg61bi2nDBrjoIrj7\n7rojkyRJS6TXyzYsnrV8NRqwzz7wve9tf4y3JEkaa17zLHVrYgLe8Q64/vq6I5EkSSPC4lnL2zvf\nafEsSZIq87INLW9bthQ3ED74ILz85XVHI0mSBszLNqRe7LQTHH88fPObdUciSZJGwMCL54h4SUSs\nj4iTB/1eUldOOw2uu67uKCRJ0ggY+GUbEXEJ8DPg3sy8sc02Xrah+jz7LBx0EDz5JOy8c93RSJKk\nAVqSyzYi4nMR8VRE3NnSflJE3BcR90fEJ+bZ73jgXuC/ga6DlAZqjz3g8MPhllvqjkSSJA25Smee\nI+IY4Hngi5n5urJtArgfOA54AlgPnJ6Z90XEmcARwK7AT4HDgF9k5rvaHN8zz6rX5ZfDPff4xEFJ\nksbckj0kJSIOAL7eVDy/Abg4M99WLp8PZGZ+ap59PwA842UbGloPPQRHHQVPPAErVtQdjSRJGpBe\ni+fJHt57X2Bz0/JjwJHzbZiZX1zsYNPT0y/OT01NMTU11UNoUocOOqgYsu4HP4Cjj647GkmS1Ccz\nMzPMzMz07Xi9nHl+N3BiZp5TLr8fODIzz+s4CM88axhcdBE88gh86EOwyy7FtNdesPvudUcmSZL6\npM4zz48D+zct71e2SaPp7LPhvPOKIvrnPy+mp5+G9evh4IPrjk6SJA2BTs48r6Y487y2XF4B/Jji\nhsH/Am4HzsjMjR0H4ZlnDavLLoMbb4RvfQvCAWMkSRp1SzVU3ZeA24A1EbEpIs7KzFngXOBm4B7g\nmm4KZ2mofexj8Pzz8NnP1h2JJEkaAgN/SEqlICLy4osv9kZBDae774a3vAU2bIBXvKLuaCRJUhfm\nbhy85JJLlmaoukHysg0NvU9+shiJ4xvf8PINSZJG2JJctiEte+efD5s3w9VX1x2JJEmqkWeeparu\nuANOOgmOOQZ22glWrYKdd4YjjoBTToF99607QkmStIgle8LgIFk8a2Rs2ACbNsGWLcX0/PNw221w\n002w//7w9rfDoYcWTymcmFh4WmybxdZH9O+1XZskSWPG4lkaBtu2wfe/X1wT/eij0Gi0n2ZnF15f\nZZvZWcgspkajv69zE2wvotsV2lXnW4/V7eRxPM6wHMvj9Odvp1SDOh+S0lfT09OOtqHRNTkJxx5b\nTOOitTjvdH7utflYvUwep57jtHaqljqeUcjRcjtOs6Xo8ER03mEf5n1GJc5u9jnuONh77+6/dwas\nX4/p9syzJEnqzFJ1nJqnTjvyS7VPHe85rPtccAGsXVvvZ7OCsTnzLEmSRoSXbmgZc6g6SZIkqSKL\nZ0mSJKkii2dJkiSpIotnSZIkqSKLZ0mSJKkii2dJkiSpoqEpnqenp/sycLUkSZLUamZmhunp6Z6P\n40NSJEmStGz0+pCUoTnzLEmSJA07i2dJkiSpIotnSZIkqSKLZ0mSJKkii2dJkiSpIotnSZIkqSKL\nZ0mSJKkii2dJkiSpIotnSZIkqSKLZ0mSJKkii2dJkiSpIotnSZIkqSKLZ0mSJKmioSmep6enmZmZ\nqTsMSZIkjaGZmRmmp6d7Pk5kZu/R9BpERA5DHJIkSRpvEUFmRrf7D82ZZ0mSJGnYWTxLkiRJFVk8\nS5IkSRVZPEuSJEkVWTxLkiRJFVk8S5IkSRVZPEuSJEkVWTxLkiRJFQ20eI6IN0fEdyLinyPiTYN8\nL0mSJGnQBn3mOYGfAauAxwb8Xir5mPP+MZf9ZT77y3z2l/nsH3PZX+ZzuFQqniPicxHxVETc2dJ+\nUkTcFxH3R8QnWvfLzO9k5inA+cAn+xOyFuN/sv4xl/1lPvvLfPaX+ewfc9lf5nO4VD3zfCVwYnND\nREwA/1C2HwacERGvLtedGRF/HxG/XW7+f8DK/oQsSZIk1WOyykaZ+d2IOKCl+Ujggcx8FCAirgFO\nA+7LzKuAqyLiXRFxIrAbRaEtSZIkjazIzGobFsXz1zPzdeXyu4ETM/Occvn9wJGZeV7HQURUC0KS\nJEnqUWZGt/tWOvM8aL38AyRJkqSl0stoG48D+zct71e2SZIkSWOpk+I5ymnOeuCVEXFARKwETgdu\n6GdwkiRJ0jCpOlTdl4DbgDURsSkizsrMWeBc4GbgHuCazNw4uFAlSZKkelUqnjPzvZm5T2auysz9\nM/PKsv3fMvNVmXlIZl7a6ZsvNk60FhYR+0XErRFxT0TcFRHnle0vi4ibI+LHEXFTROxWd6yjIiIm\nImJDRNxQLpvLLkXEbhFxbURsLD+jv2c+uxcRfxYRd0fEnRFxdUSsNJ/Vzfe8goXyFxEXRMQD5ef3\n9+uJeni1yefflvn6UUR8LSJ2bVpnPtto9yyNct3HI6IREbs3tZnLBSzwbJJzy5zdFRGXNrV3nM9B\nP2GwrYXGiVZl24A/z8zDgKOAPylzeD5wS2a+CrgVuKDGGEfNR4F7m5bNZfc+DdyYma8Bfge4D/PZ\nlYjYh+KXviPKEY8mgTMwn534tecV0CZ/EXEo8B7gNcDbgH+KCG9s39F8+bwZOCwzXw88gPmsar5c\nEhH7AScAjza1vQZzuZj5nk0yBbwDWJuZa4HLyvau8llb8UzTONGZ+QIwN060KsrMJzPzR+X888BG\nihs3TwO+UG72BeCd9UQ4Wso/VCcDVzQ1m8sulGecjm36lWpbZv4U89mLFcAuETEJ7Exxg7b5rCgz\nvwv8b0tzu/ydSnEp4rbMfISiEDxyKeIcFfPlMzNvycxGubiO4vsIzOeC2nw2AS4H/qKl7TTM5YLa\n5POPgUszc1u5zTNle1f5rLN43hfY3LT8WNmmLkTEauD1FH+w9srMp6AosIE964tspMz9oWoed9xc\ndudA4JmIuLK8DOYzEfESzGdXMvMJ4O+ATRRF808z8xbMZ6/2bJO/1u+nx/H7qVNnAzeW8+azQxFx\nKrA5M+9qWWUuu7MGeFNErIuI/4iI3y3bu8pnncWz+iQiXgp8FfhoeQa69aEzPoRmERFxCvBUeSZ/\noZ9szGU1k8ARwD9m5hHAzyl+Ivez2YWI+C2KMyQHAPtQnIF+H+az38xfH0TEXwEvZOaX645lFEXE\nzsCFwMV1xzJGJoGXZeYbgL8Eru3lYHUWz44T3QflT7hfBa7KzOvL5qciYq9y/d7A03XFN0LeCJwa\nEQ8BXwbeGhFXAU+ay648RnHW5I5y+WsUxbSfze4cDzyUmf9TjnT0r8DRmM9etcvf48Armrbz+6mi\niPggxeVv721qNp+dORhYDfxnRDxMka8NEbEn1k7d2gz8C0BmrgdmI2IPusxnncWz40T3x+eBezPz\n001tNwAfLOf/CLi+dSftKDMvLEeSOYjis3hrZp4JfB1z2bHyp/DNEbGmbDqOYkhLP5vd2QS8ISJ2\nKm9mOY7ixlbz2ZnW5xW0y98NwOnliCYHAq8Ebl+qIEfIDvmMiJMoLn07NTO3Nm1nPhf3Yi4z8+7M\n3DszD8rMAylORhyemU9T5PIPzeWiWv+vXwe8FaD8XlqZmc/SZT5rezx3Zs5GxJ9S3J07AXzOcaI7\nExFvBN4H3BURP6T4yfFC4FPAVyLibIq7dN9TX5Qj71LMZbfOA66OiN8AHgLOorjpzXx2KDNvj4iv\nAj8EXihfPwP8JuazkiieVzAF7BERmyh+Er8UuLY1f5l5b0R8haKD8gLwkcz0ko4mbfJ5IbAS+Pdy\nwIJ1mfkR87mw+XI5d7N1KdleWJvLRbT5bH4euDIi7gK2Ah+A7vMZ5lySJEmqxhsGJUmSpIosniVJ\nkqSKLJ4lSZKkiiyeJUmSpIosniVJkqSKLJ4lSZKkiiyeJUmSpIr+H/cm0Q4YvCosAAAAAElFTkSu\nQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe6e5af5210>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "semilogy(sigmaw)\n",
    "semilogy(sigmau)\n",
    "semilogy(sigmav)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Use this code to identify the PDE from reduced basis\n",
    "# dim = 50\n",
    "# Wn = uw[:,0:dim].dot(np.diag(sigmaw[0:dim]).dot(vw[:,0:dim].T)).reshape(n,m,steps)\n",
    "# Un = uu[:,0:dim].dot(np.diag(sigmau[0:dim]).dot(vu[:,0:dim].T)).reshape(n,m,steps)\n",
    "# Vn = uv[:,0:dim].dot(np.diag(sigmav[0:dim]).dot(vv[:,0:dim].T)).reshape(n,m,steps)\n",
    "\n",
    "# Or this code to take the full solution\n",
    "Wn = W.reshape(n,m,steps)\n",
    "Un = U.reshape(n,m,steps)\n",
    "Vn = V.reshape(n,m,steps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sample data points\n",
    "\n",
    "Now randomly sample some points to use.  See figure 1, panel 2a-2c for visual description of what we're doing here.  5000 spatial points are samples and 60 timepoints are viewed at each one.  Note that we still need nearby points to take derivatives."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Sample a collection of data points, stay away from edges so I can just use centered finite differences.\n",
    "np.random.seed(0)\n",
    "\n",
    "num_xy = 5000\n",
    "num_t = 60\n",
    "num_points = num_xy * num_t\n",
    "boundary = 5\n",
    "boundary_x = 10\n",
    "points = {}\n",
    "count = 0\n",
    "\n",
    "for p in range(num_xy):\n",
    "    x = np.random.choice(np.arange(boundary_x,n-boundary_x),1)[0]\n",
    "    y = np.random.choice(np.arange(boundary,m-boundary),1)[0]\n",
    "    for t in range(num_t):\n",
    "        points[count] = [x,y,2*t+12]\n",
    "        count = count + 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Construct $\\Theta (U)$ and compute $U_t$\n",
    "\n",
    "Take derivatives and assemble into $\\Theta(\\omega, u ,v)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Take up to second order derivatives.\n",
    "w = np.zeros((num_points,1))\n",
    "u = np.zeros((num_points,1))\n",
    "v = np.zeros((num_points,1))\n",
    "wt = np.zeros((num_points,1))\n",
    "wx = np.zeros((num_points,1))\n",
    "wy = np.zeros((num_points,1))\n",
    "wxx = np.zeros((num_points,1))\n",
    "wxy = np.zeros((num_points,1))\n",
    "wyy = np.zeros((num_points,1))\n",
    "\n",
    "N = 2*boundary-1  # odd number of points to use in fitting\n",
    "Nx = 2*boundary_x-1  # odd number of points to use in fitting\n",
    "deg = 5 # degree of polynomial to use\n",
    "\n",
    "for p in points.keys():\n",
    "    \n",
    "    [x,y,t] = points[p]\n",
    "    w[p] = Wn[x,y,t]\n",
    "    u[p] = Un[x,y,t]\n",
    "    v[p] = Vn[x,y,t]\n",
    "    \n",
    "    wt[p] = PolyDiffPoint(Wn[x,y,t-(N-1)/2:t+(N+1)/2], np.arange(N)*dt, deg, 1)[0]\n",
    "    \n",
    "    x_diff = PolyDiffPoint(Wn[x-(Nx-1)/2:x+(Nx+1)/2,y,t], np.arange(Nx)*dx, deg, 2)\n",
    "    y_diff = PolyDiffPoint(Wn[x,y-(N-1)/2:y+(N+1)/2,t], np.arange(N)*dy, deg, 2)\n",
    "    wx[p] = x_diff[0]\n",
    "    wy[p] = y_diff[0]\n",
    "    \n",
    "    x_diff_yp = PolyDiffPoint(Wn[x-(Nx-1)/2:x+(Nx+1)/2,y+1,t], np.arange(Nx)*dx, deg, 2)\n",
    "    x_diff_ym = PolyDiffPoint(Wn[x-(Nx-1)/2:x+(Nx+1)/2,y-1,t], np.arange(Nx)*dx, deg, 2)\n",
    "    \n",
    "    wxx[p] = x_diff[1]\n",
    "    wxy[p] = (x_diff_yp[0]-x_diff_ym[0])/(2*dy)\n",
    "    wyy[p] = y_diff[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Candidate terms for PDE\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "['1',\n",
       " 'w_{x}',\n",
       " 'w_{y}',\n",
       " 'w_{xx}',\n",
       " 'w_{xy}',\n",
       " 'w_{yy}',\n",
       " 'uv',\n",
       " 'v^2',\n",
       " 'u^2',\n",
       " 'w',\n",
       " 'wu',\n",
       " 'v',\n",
       " 'wv',\n",
       " 'u',\n",
       " 'w^2',\n",
       " 'uvw_{x}',\n",
       " 'v^2w_{x}',\n",
       " 'u^2w_{x}',\n",
       " 'ww_{x}',\n",
       " 'wuw_{x}',\n",
       " 'vw_{x}',\n",
       " 'wvw_{x}',\n",
       " 'uw_{x}',\n",
       " 'w^2w_{x}',\n",
       " 'uvw_{y}',\n",
       " 'v^2w_{y}',\n",
       " 'u^2w_{y}',\n",
       " 'ww_{y}',\n",
       " 'wuw_{y}',\n",
       " 'vw_{y}',\n",
       " 'wvw_{y}',\n",
       " 'uw_{y}',\n",
       " 'w^2w_{y}',\n",
       " 'uvw_{xx}',\n",
       " 'v^2w_{xx}',\n",
       " 'u^2w_{xx}',\n",
       " 'ww_{xx}',\n",
       " 'wuw_{xx}',\n",
       " 'vw_{xx}',\n",
       " 'wvw_{xx}',\n",
       " 'uw_{xx}',\n",
       " 'w^2w_{xx}',\n",
       " 'uvw_{xy}',\n",
       " 'v^2w_{xy}',\n",
       " 'u^2w_{xy}',\n",
       " 'ww_{xy}',\n",
       " 'wuw_{xy}',\n",
       " 'vw_{xy}',\n",
       " 'wvw_{xy}',\n",
       " 'uw_{xy}',\n",
       " 'w^2w_{xy}',\n",
       " 'uvw_{yy}',\n",
       " 'v^2w_{yy}',\n",
       " 'u^2w_{yy}',\n",
       " 'ww_{yy}',\n",
       " 'wuw_{yy}',\n",
       " 'vw_{yy}',\n",
       " 'wvw_{yy}',\n",
       " 'uw_{yy}',\n",
       " 'w^2w_{yy}']"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Form a huge matrix using up to quadratic polynomials in all variables.\n",
    "X_data = np.hstack([w,u,v])\n",
    "X_ders = np.hstack([np.ones((num_points,1)), wx, wy, wxx, wxy, wyy])\n",
    "X_ders_descr = ['','w_{x}', 'w_{y}','w_{xx}','w_{xy}','w_{yy}']\n",
    "X, description = build_Theta(X_data, X_ders, X_ders_descr, 2, data_description = ['w','u','v'])\n",
    "print 'Candidate terms for PDE'\n",
    "['1']+description[1:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve for $\\xi$\n",
    "\n",
    "TrainSTRidge splits the data up into 80% for training and 20% for validation.  It searches over various tolerances in the STRidge algorithm and finds the one with the best performance on the validation set, including an $\\ell^0$ penalty for $\\xi$ in the loss function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "w_t = (0.009884 +0.000000i)w_{xx}\n",
      "    + (0.009902 +0.000000i)w_{yy}\n",
      "    + (-0.990371 +0.000000i)uw_{x}\n",
      "    + (-0.986629 +0.000000i)vw_{y}\n",
      "   \n"
     ]
    }
   ],
   "source": [
    "lam = 10**-5\n",
    "d_tol = 5\n",
    "c = TrainSTRidge(X,wt,lam,d_tol)\n",
    "print_pde(c, description, ut = 'w_t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error using PDE-FIND to identify Navier-Stokes:\n",
      "\n",
      "Mean parameter error: 1.11 %\n",
      "Standard deviation of parameter error: 0.152161772466 %\n"
     ]
    }
   ],
   "source": [
    "err = abs(np.array([(0.009884-0.01)*100/0.01, (0.009902-0.01)*100/0.01, (-0.990371+1)*100, (-0.986629+1)*100]))\n",
    "print \"Error using PDE-FIND to identify Navier-Stokes:\\n\"\n",
    "print \"Mean parameter error:\", np.mean(err), '%'\n",
    "print \"Standard deviation of parameter error:\", np.std(err), '%'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
